/*********************************************************
 *
 * MDCM
 *
 * Modifications github.com/issakomi
 *
 *********************************************************/

/*=========================================================================

  Program: GDCM (Grassroots DICOM). A DICOM library

  Copyright (c) 2006-2011 Mathieu Malaterre
  All rights reserved.
  See Copyright.txt or http://gdcm.sourceforge.net/Copyright.html for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

#include "mdcmIPPSorter.h"
#include "mdcmScanner.h"
#include "mdcmElement.h"
#include "mdcmDirectionCosines.h"
#include <map>
#include <cmath>

namespace mdcm
{

IPPSorter::IPPSorter()
{
  ComputeZSpacing = true;
  DropDuplicatePositions = false;
  ZSpacing = 0;
  ZTolerance = 1e-6;
  DirCosTolerance = 0.;
}


inline double spacing_round(double n, int d) /* pow is defined as pow(double, double) or pow(double int) on M$ comp */
{
  return floor(n * pow(10., d) + .5) / pow(10., d);
}

bool IPPSorter::Sort(std::vector<std::string> const & filenames)
{
  // BUG: I cannot clear Filenames since input filenames could also be the output of ourself
  ZSpacing = 0;
  if(filenames.empty())
  {
    Filenames.clear();
    return true;
  }
  Scanner scanner;
  const Tag tipp(0x0020,0x0032);    // Image Position (Patient)
  const Tag tiop(0x0020,0x0037);    // Image Orientation (Patient)
  const Tag tframe(0x0020,0x0052);  // Frame of Reference UID
  const Tag tgantry(0x0018,0x1120); // Gantry/Detector Tilt
  // Temporal Position Identifier (0020,0100) 3 Temporal order of a dynamic or functional set of Images.
  //const Tag tpi(0x0020,0x0100);
  scanner.AddTag(tipp);
  scanner.AddTag(tiop);
  scanner.AddTag(tframe);
  scanner.AddTag(tgantry);
  bool b = scanner.Scan(filenames);
  if(!b)
  {
    mdcmDebugMacro("Scanner failed");
    return false;
  }
  Scanner::ValuesType gantry = scanner.GetValues(tgantry);
  if(gantry.size() > 1)
  {
    mdcmDebugMacro("More than one Gantry/Detector Tilt");
    return false;
  }
  if(gantry.size() == 1)
  {
    std::stringstream ss(*gantry.begin());
    double tilt;
    ss >> tilt;
    if(tilt != 0.0)
    {
      mdcmDebugMacro("Gantry/Detector Tilt is not 0");
      return false;
    }
  }
  Scanner::ValuesType iops = scanner.GetValues(tiop);
  Scanner::ValuesType frames = scanner.GetValues(tframe);
  if(DirCosTolerance == 0.)
  {
    if(iops.size() != 1)
    {
      mdcmDebugMacro("More than one IOP (or no IOP): " << iops.size());
      //std::copy(iops.begin(), iops.end(), std::ostream_iterator<std::string>(std::cout, "\n"));
      return false;
    }
  }
  const size_t fsize = frames.size(); // Should I really tolerate issue with Frame of Reference UID ?
  if(fsize == 1) // by the book
  {
    // TODO: need to check not empty ? technically PMS used to send MR Image Storage with empty FoR
  }
  else if(fsize == 0 || fsize == filenames.size()) // Should I really tolerate no Frame of Reference UID ?
  {
    mdcmWarningMacro("Odd number of Frame Of Reference UID (continuing with caution): " << fsize);
  }
  else
  {
    mdcmErrorMacro("Sorry your setup with Frame Of Reference UID does not make any sense: " << fsize);
    return false;
  }
  const char *reference = filenames[0].c_str();
  // we cannot simply consider the first file, what if this is not DICOM ?
  for(std::vector<std::string>::const_iterator it1 = filenames.begin();
    it1 != filenames.end(); ++it1)
  {
    const char *filename = it1->c_str();
    bool iskey = scanner.IsKey(filename);
    if(iskey)
    {
      reference = filename;
    }
  }
  Scanner::TagToValue const &t2v = scanner.GetMapping(reference);
  Scanner::TagToValue::const_iterator it = t2v.find(tiop);
  // Take the first file in the list of filenames, if not IOP is found, simply gives up:
  if(it == t2v.end())
  {
    // Dead code
    mdcmDebugMacro("No iop in: " << reference);
    return false;
  }
  if(it->first != tiop)
  {
    // first file does not contains Image Orientation (Patient), let's give up
    mdcmDebugMacro("No iop in first file ");
    return false;
  }
  const char *dircos = it->second;
  if(!dircos)
  {
    // first file does contains Image Orientation (Patient), but it is empty
    mdcmDebugMacro("Empty iop in first file ");
    return false;
  }
  // http://www.itk.org/pipermail/insight-users/2003-September/004762.html
  // Compute normal:
  // The steps I take when reconstructing a volume are these: First,
  // calculate the slice normal from IOP:
  double normal[3];
  DirectionCosines dc;
  dc.SetFromString(dircos);
  if(!dc.IsValid()) return false;
  dc.Cross(normal);
  // You only have to do this once for all slices in the volume. Next, for
  // each slice, calculate the distance along the slice normal using the IPP
  // tag ("dist" is initialized to zero before reading the first slice) :
  //typedef std::multimap<double, const char*> SortedFilenames;
  typedef std::map<double, const char*> SortedFilenames;
  SortedFilenames sorted;
  {
    std::vector<std::string>::const_iterator it1 = filenames.begin();
    DirectionCosines dc2;
    for(; it1 != filenames.end(); ++it1)
    {
      const char *filename = it1->c_str();
      bool iskey = scanner.IsKey(filename);
      if(iskey)
      {
        const char *value =  scanner.GetValue(filename, tipp);
        if(value)
        {
          if(DirCosTolerance != 0.)
          {
            const char *value2 =  scanner.GetValue(filename, tiop);
            if(!dc2.SetFromString(value2))
            {
              if(value2)
                mdcmWarningMacro(filename << " cant read IOP: " << value2);
              return false;
            }
            double cd = dc2.CrossDot(dc);
            // result should be as close to 1 as possible:
            if(fabs(1 - cd) > DirCosTolerance)
            {
              mdcmWarningMacro(filename << " Problem with DirCosTolerance: ");
              // Cant print cd since 0.9999 is printed as 1... may confuse user
              return false;
            }
          }
          //mdcmDebugMacro(filename << " has " << ipp << " = " << value);
          Element<VR::DS,VM::VM3> ipp;
          std::stringstream ss;
          ss.str(value);
          ipp.Read(ss);
          double dist = 0;
          for (int i = 0; i < 3; ++i) dist += normal[i]*ipp[i];
          // FIXME: This test is weak, since implicitely we are doing a != on floating point value
          if(sorted.find(dist) != sorted.end())
          {
            if(this->DropDuplicatePositions)
            {
              mdcmWarningMacro("dropping file " << filename << " since Z position: " << dist << " already found");
              continue;
            }
            mdcmWarningMacro("dist: " << dist << " for " << filename <<
              " already found in " << sorted.find(dist)->second);
            return false;
          }
          sorted.insert(SortedFilenames::value_type(dist,filename));
        }
        else
        {
          mdcmDebugMacro("File: " << filename << " has no Tag" << tipp << ". Skipping.");
        }
      }
      else
      {
        mdcmDebugMacro("File: " << filename << " could not be read. Skipping.");
      }
    }
  }
  assert(!sorted.empty());
  {
    SortedFilenames::const_iterator it2 = sorted.begin();
    double prev = it2->first;
    Filenames.push_back(it2->second);
    if(sorted.size() > 1)
    {
      bool spacingisgood = true;
      ++it2;
      double current = it2->first;
      double zspacing = current - prev;
      for(; it2 != sorted.end(); ++it2)
      {
        //std::cout << it2->first << " " << it2->second << std::endl;
        current = it2->first;
        Filenames.push_back(it2->second);
        if(fabs((current - prev) - zspacing) > ZTolerance)
        {
          mdcmDebugMacro("ZTolerance test failed. You need to decrease ZTolerance.");
          spacingisgood = false;
        }
        // update prev for the next for-loop
        prev = current;
      }
      // is spacing good ?
      if(spacingisgood && ComputeZSpacing)
      {
        // If user ask for a ZTolerance of 1e-4, there is no need for us to
        // store the extra digits... this will make sure to return 2.2 from a 2.1999938551239993 value
        const int l = (int)(-log10(ZTolerance));
        ZSpacing = spacing_round(zspacing, l);
      }
      if(!spacingisgood)
      {
        std::ostringstream os;
        os << "Filenames and 'z' positions" << std::endl;
        double prev1 = 0.;
        for(SortedFilenames::const_iterator it1 = sorted.begin(); it1 != sorted.end(); ++it1)
        {
          std::string f = it1->second;
          if(f.length() > 62)
          {
            f = f.substr(0,10) + " ... " + f.substr(f.length()-47);
          }
          double d = it1->first - prev1;
          if(it1 != sorted.begin() && fabs(d - zspacing) > ZTolerance) os << "* ";
          else os << "  ";
          os << it1->first << "\t" << f << std::endl;
          prev1 = it1->first;
        }
        mdcmDebugMacro(os.str());
      }
      assert(spacingisgood == false ||  (ComputeZSpacing ? (ZSpacing > ZTolerance && ZTolerance > 0) : ZTolerance > 0));
    }
  }
  // return true: means sorting succeed, it does not mean spacing computation succeded
  return true;
}

} // end namespace mdcm
